#include "Process.h"
#include "Blob.h"

#include "EdgesSubPix.h"

//#include "frangi.h"

using namespace std;
using namespace cv;

stop_watch  tc;


#include <iostream>
#include <vector>
#include <sstream>
#include <fstream>


//#include "lshaped_fitting.h"

int main(int argc, char** argv)
{

    //// 3D Cloud --> 2D Dimension

    //// Every Cluster Points.
    //std::vector<cv::Point2f> hull;

    //// Load Point Cloud For Shaped-BBox Fit.
    //hull.push_back(cv::Point2f(19.61026979843742, 7.52767729296133));
    //hull.push_back(cv::Point2f(17.03273076303270, 7.58346032483931));
    //hull.push_back(cv::Point2f(15.13402863573614, 7.71117274328501));
    //hull.push_back(cv::Point2f(14.72757894312114, 8.50297166732246));
    //hull.push_back(cv::Point2f(14.77044667030188, 9.59204022261293));
    //hull.push_back(cv::Point2f(14.76137992714666, 10.7247702891167));
    //hull.push_back(cv::Point2f(14.74829228186517, 11.9429316067475));

    //// Do Shaped-BBox Fit.
    //LShapedFIT lshaped;

    //cv::RotatedRect rr = lshaped.FitBox(&hull);
    //std::cout << "Shaped-BBox Message : " << rr.size.width << " " << rr.size.height << " " << rr.angle;

    //// --- Vertex Standard Output Coordinates.
    //// [19.60907603446491,  7.493492365484474]
    //// [19.75834431780993,  11.76797673445868]
    //// [14.74829228186517,  11.94293160674757]
    //// [14.59902399852014,  7.668447237773362]
    //std::vector<cv::Point2f> vertices = lshaped.getRectVertex();

    //std::cout << "Top 4 Vertices" << std::endl;
    //for (size_t i = 0; i < vertices.size(); ++i)
    //    std::cout << "  " << vertices[i].x << "  " << vertices[i].y << std::endl;

    //return 0;
}


template <typename T>
int sgn(T val) {

    T tt = (T(0) < val);
    T ttm = (val < T(0));
    return (T(0) < val) - (val < T(0));
}
void exportMesh(cv::Mat Depth, cv::Mat Normals, cv::Mat texture) {
    /* writing obj for export */
    std::ofstream objFile, mtlFile;
    objFile.open("export.obj");
    int width = Depth.cols;
    int height = Depth.rows;
    /* vertices, normals, texture coords */
    objFile << "mtllib export.mtl" << std::endl;

    for (int y = 0; y < height; y++) {
        for (int x = 0; x < width; x++) {
            objFile << "v " << x << " " << y << " "
                << Depth.at<float>(cv::Point(x, y)) << std::endl;
            objFile << "vt " << x / (width - 1.0f) << " " << (1.0f - y) / height
                << " "
                << "0.0" << std::endl;
            objFile << "vn " << (float)Normals.at<cv::Vec3b>(y, x)[0] << " "
                << (float)Normals.at<cv::Vec3b>(y, x)[1] << " "
                << (float)Normals.at<cv::Vec3b>(y, x)[2] << std::endl;
        }
    }

    /* faces */
    objFile << "usemtl picture" << std::endl;

    for (int y = 0; y < height - 1; y++) {
        for (int x = 0; x < width - 1; x++) {
            int f1 = x + y * width + 1;
            int f2 = x + y * width + 2;
            int f3 = x + (y + 1) * width + 1;
            int f4 = x + (y + 1) * width + 2;
            objFile << "f " << f1 << "/" << f1 << "/" << f1 << " ";
            objFile << f2 << "/" << f2 << "/" << f2 << " ";
            objFile << f3 << "/" << f3 << "/" << f3 << std::endl;
            objFile << "f " << f2 << "/" << f2 << "/" << f2 << " ";
            objFile << f4 << "/" << f4 << "/" << f4 << " ";
            objFile << f3 << "/" << f3 << "/" << f3 << std::endl;
        }
    }

    /* texture */
    cv::imwrite("export.jpg", texture);
    mtlFile.open("export.mtl");
    mtlFile << "newmtl picture" << std::endl;
    mtlFile << "map_Kd export.jpg" << std::endl;
    objFile.close();
    mtlFile.close();
}

cv::Mat imageMask(std::vector<cv::Mat> camImages) 
{
    assert(camImages.size() > 0);
    cv::Mat image = camImages[0].clone();
    int quarter = image.cols / 4.0;
    int eighth = image.rows / 8.0;
    cv::Mat result, bgModel, fgModel;
    cv::Rect area(quarter, eighth, 3 * quarter, 7 * eighth);
    /* grabcut expects rgb images */
    cv::cvtColor(image, image, COLOR_GRAY2BGR);
    cv::grabCut(image, result, area, bgModel, fgModel, 1, cv::GC_INIT_WITH_RECT);
    cv::compare(result, cv::GC_PR_FGD, result, cv::CMP_EQ);
    return result;
}

/**
 * Returns (binary) light pattern img L_j {j = 1..N}
 */
cv::Mat lightPattern(int width, int height, int j, int N) 
{
    cv::Mat img(height, width, CV_8UC1, cv::Scalar::all(0));

    for (int y = -(height / 2); y < height / 2; y++) {
        for (int x = -(width / 2); x < width / 2; x++) {
            if (sgn(x * cos(2 * CV_PI * j / N) +
                y * sin(2 * CV_PI * j / N)) == 1) {
                img.at<uchar>(y + height / 2, x + width / 2) = 255;
            }
        }
    }

    return img;
}

cv::Mat computeNormals(std::vector<cv::Mat> camImages,cv::Mat Mask = cv::Mat()) 
{
    int height = camImages[0].rows;
    int width = camImages[0].cols;
    int numImgs = camImages.size();
    /* populate A */
    cv::Mat A(height * width, numImgs, CV_32FC1, cv::Scalar::all(0));
    for (int k = 0; k < numImgs; k++) 
    {
        int idx = 0;
        for (int i = 0; i < height; i++) 
        {
            for (int j = 0; j < width; j++) 
            {
                A.at<float>(idx++, k) = camImages[k].data[i * width + j];
            }
        }
    }

    /* speeding up computation, SVD from A^TA instead of AA^T */
    cv::Mat U, S, Vt;
    cv::SVD::compute(A.t(), S, U, Vt, cv::SVD::MODIFY_A);
    cv::Mat EV = Vt.t();
    cv::Mat N(height, width, CV_8UC3, cv::Scalar::all(0));
    int idx = 0;

    for (int i = 0; i < height; i++) 
    {
        for (int j = 0; j < width; j++) 
        {
            float rSxyz = 1.0f / sqrt(EV.at<float>(idx, 0) * EV.at<float>(idx, 0) +
                EV.at<float>(idx, 1) * EV.at<float>(idx, 1) +
                EV.at<float>(idx, 2) * EV.at<float>(idx, 2));
            /* V contains the eigenvectors of A^TA, which are as well the z,x,y
             * components of the surface normals for each pixel	*/
            float sz = 128.0f +
                127.0f * sgn(EV.at<float>(idx, 0)) *
                fabs(EV.at<float>(idx, 0)) * rSxyz;
            float sx = 128.0f +
                127.0f * sgn(EV.at<float>(idx, 1)) *
                fabs(EV.at<float>(idx, 1)) * rSxyz;
            float sy = 128.0f +
                127.0f * sgn(EV.at<float>(idx, 2)) *
                fabs(EV.at<float>(idx, 2)) * rSxyz;
            N.at<cv::Vec3b>(i, j) = cv::Vec3b(sx, sy, sz);

            idx += 1;
        }
    }

    return N;
}


void NormalsSpeed(std::vector<cv::Mat> camImages, cv::Mat& XVal, cv::Mat& YVal, cv::Mat& ZVal)
{
    int height = camImages[0].rows;
    int width = camImages[0].cols;
    int numImgs = camImages.size();
    cv::Mat A( numImgs, height * width, CV_32FC1);

    static int AS = 16;//AVX2
    int LastStart = width - AS;
    size_t End = AlignLo(width, AS);//AVX2
    bool EndType = End != width;

    uchar* SrcPtr;
    float* DisPtr;
    float* DisPtrStart;



    for (int k = 0; k < numImgs; k++)
    {
        DisPtr = A.ptr<float>(k);
        for (int i = 0; i < height; i++)
        {
            SrcPtr = camImages[k].ptr<uchar>(i);         
            for (int j = 0; j < End; j += AS)
            {
                __m128i SrcD = _mm_loadu_epi8(SrcPtr + j);
                _mm256_storeu_ps(DisPtr + j, _mm256_cvtepi32_ps(_mm256_cvtepu8_epi32(SrcD)));
                _mm256_storeu_ps(DisPtr + 8 + j, _mm256_cvtepi32_ps(_mm256_cvtepu8_epi32(_mm_srli_si128(SrcD, 8))));
            }
            if (EndType)
            {
                __m128i SrcD = _mm_loadu_epi8(SrcPtr + LastStart);
                _mm256_storeu_ps(DisPtr + LastStart, _mm256_cvtepi32_ps(_mm256_cvtepu8_epi32(SrcD)));
                _mm256_storeu_ps(DisPtr + LastStart + 8, _mm256_cvtepi32_ps(_mm256_cvtepu8_epi32(_mm_srli_si128(SrcD, 8))));
            }
            DisPtr += width;
        }

    }


    cv::Mat U, S, Vt;
    cv::SVD::compute(A, S, U, Vt, cv::SVD::MODIFY_A);
    XVal = cv::Mat(height, width, CV_8UC1);
    YVal = cv::Mat(height, width, CV_8UC1);
    ZVal = cv::Mat(height, width, CV_8UC1);
    AS = 16;//AVX2
    LastStart = width - AS;
    End = AlignLo(width, AS);//AVX2
    EndType = End != width;


    uchar* XPtr;
    uchar* YPtr;
    uchar* ZPtr;
    float* SrcPtrA = Vt.ptr<float>(0);
    float* SrcPtrB = Vt.ptr<float>(1);
    float* SrcPtrC = Vt.ptr<float>(2);
    __m256 valA, valB, valC;
    __m256 OnePs = _mm256_set1_ps(1);
    __m256 Stac128 = _mm256_set1_ps(128);
    __m256 Stac127 = _mm256_set1_ps(127);
    __m256 szh, syh, sxh;
    __m256 szl, syl, sxl;
    for (int i = 0; i < height; i++)
    {   
        XPtr = XVal.ptr<uchar>(i);
        YPtr = YVal.ptr<uchar>(i);
        ZPtr = ZVal.ptr<uchar>(i);
        for (int j = 0; j < End; j+= AS)
        {
            {
                valA = _mm256_loadu_ps(SrcPtrA + j);
                valB = _mm256_loadu_ps(SrcPtrB + j);
                valC = _mm256_loadu_ps(SrcPtrC + j);
                __m256 rSxyz = _mm256_div_ps(OnePs, _mm256_sqrt_ps(_mm256_add_ps(_mm256_add_ps(_mm256_mul_ps(valA, valA), _mm256_mul_ps(valB, valB)), _mm256_mul_ps(valC, valC))));
                szl = _mm256_add_ps(Stac128, _mm256_mul_ps(_mm256_mul_ps(Stac127, valA), rSxyz));
                sxl = _mm256_add_ps(Stac128, _mm256_mul_ps(_mm256_mul_ps(Stac127, valB), rSxyz));
                syl = _mm256_add_ps(Stac128, _mm256_mul_ps(_mm256_mul_ps(Stac127, valC), rSxyz));
            }
            {
                valA = _mm256_loadu_ps(SrcPtrA + j + 8);
                valB = _mm256_loadu_ps(SrcPtrB + j + 8);
                valC = _mm256_loadu_ps(SrcPtrC + j + 8);
                __m256 rSxyz = _mm256_div_ps(OnePs, _mm256_sqrt_ps(_mm256_add_ps(_mm256_add_ps(_mm256_mul_ps(valA, valA), _mm256_mul_ps(valB, valB)), _mm256_mul_ps(valC, valC))));
                szh = _mm256_add_ps(Stac128, _mm256_mul_ps(_mm256_mul_ps(Stac127, valA), rSxyz));
                sxh = _mm256_add_ps(Stac128, _mm256_mul_ps(_mm256_mul_ps(Stac127, valB), rSxyz));
                syh = _mm256_add_ps(Stac128, _mm256_mul_ps(_mm256_mul_ps(Stac127, valC), rSxyz));
            }
            _mm_storeu_epi8(ZPtr + j, CvtPsToEpi8(szh, szl));
            _mm_storeu_epi8(XPtr + j, CvtPsToEpi8(sxh, sxl));
            _mm_storeu_epi8(YPtr + j, CvtPsToEpi8(syh, syl));
        }
        if (EndType)
        {
            {
                valA = _mm256_loadu_ps(SrcPtrA + LastStart);
                valB = _mm256_loadu_ps(SrcPtrB + LastStart);
                valC = _mm256_loadu_ps(SrcPtrC + LastStart);
                __m256 rSxyz = _mm256_div_ps(OnePs, _mm256_sqrt_ps(_mm256_add_ps(_mm256_add_ps(_mm256_mul_ps(valA, valA), _mm256_mul_ps(valB, valB)), _mm256_mul_ps(valC, valC))));
                szl = _mm256_add_ps(Stac128, _mm256_mul_ps(_mm256_mul_ps(Stac127, valA), rSxyz));
                sxl = _mm256_add_ps(Stac128, _mm256_mul_ps(_mm256_mul_ps(Stac127, valB), rSxyz));
                syl = _mm256_add_ps(Stac128, _mm256_mul_ps(_mm256_mul_ps(Stac127, valC), rSxyz));
            }
            {
                valA = _mm256_loadu_ps(SrcPtrA + LastStart + 8);
                valB = _mm256_loadu_ps(SrcPtrB + LastStart + 8);
                valC = _mm256_loadu_ps(SrcPtrC + LastStart + 8);
                __m256 rSxyz = _mm256_div_ps(OnePs, _mm256_sqrt_ps(_mm256_add_ps(_mm256_add_ps(_mm256_mul_ps(valA, valA), _mm256_mul_ps(valB, valB)), _mm256_mul_ps(valC, valC))));
                szh = _mm256_add_ps(Stac128, _mm256_mul_ps(_mm256_mul_ps(Stac127, valA), rSxyz));
                sxh = _mm256_add_ps(Stac128, _mm256_mul_ps(_mm256_mul_ps(Stac127, valB), rSxyz));
                syh = _mm256_add_ps(Stac128, _mm256_mul_ps(_mm256_mul_ps(Stac127, valC), rSxyz));
            }
            _mm_storeu_epi8(ZPtr + LastStart, CvtPsToEpi8(szh, szl));
            _mm_storeu_epi8(XPtr + LastStart, CvtPsToEpi8(sxh, sxl));
            _mm_storeu_epi8(YPtr + LastStart, CvtPsToEpi8(syh, syl));
        }
        SrcPtrA += width;
        SrcPtrB += width;
        SrcPtrC += width;
    }
}

void HeightsSpeed(cv::Mat& XVal, cv::Mat& YVal, cv::Mat& Z, int iterations)
{
    static int AS = 16;//AVX2
    int LastStart = XVal.cols - 1 - AS;
    size_t End = AlignLo(XVal.cols - 2, AS);//AVX2
    bool EndType = End != XVal.cols - 2;
    End += 1;

    uchar* XPtrU;
    uchar* YPtrU;
    uchar* XPtrD;
    uchar* YPtrD;
    uchar* XPtr;
    uchar* YPtr;

    float* SrcU;
    float* SrcD;
    float* Src;

    __m128i Zero128 = _mm_set1_epi8(0);
    __m256 val025 = _mm256_set1_ps(0.25);
    __m256 valZero = _mm256_set1_ps(0);

    for (int k = 0; k < iterations; k++)
    {
        for (int i = 1; i < XVal.rows - 1; i++)
        {
            XPtrU = XVal.ptr<uchar>(i - 1);
            YPtrU = YVal.ptr<uchar>(i - 1);
            XPtrD = XVal.ptr<uchar>(i + 1);
            YPtrD = YVal.ptr<uchar>(i + 1);
            XPtr = XVal.ptr<uchar>(i);
            YPtr = YVal.ptr<uchar>(i);

            SrcU = Z.ptr<float>(i - 1);
            SrcD = Z.ptr<float>(i + 1);
            Src = Z.ptr<float>(i);
            for (int j = 1; j < End; j++)
            {
                __m128i valXU = _mm_loadu_epi8(XPtrU + j);//
                __m128i valYL = _mm_loadu_epi8(YPtr + j - 1);//
                __m128i cenX = _mm_loadu_epi8(XPtr + j);//
                __m128i cenY = _mm_loadu_epi8(YPtr + j);//

                __m256 valA = _mm256_add_ps(_mm256_add_ps(_mm256_add_ps(_mm256_loadu_ps(SrcU + j), _mm256_loadu_ps(SrcD + j)), _mm256_loadu_ps(Src + j - 1)), _mm256_loadu_ps(Src + j + 1));
                __m256 valB = _mm256_sub_ps(_mm256_cvtepi32_ps(_mm256_cvtepu8_epi32(valXU)), _mm256_cvtepi32_ps(_mm256_cvtepu8_epi32(cenX)));
                __m256 valC = _mm256_sub_ps(_mm256_cvtepi32_ps(_mm256_cvtepu8_epi32(valYL)), _mm256_cvtepi32_ps(_mm256_cvtepu8_epi32(cenY)));
                __m256 vallo = _mm256_mul_ps(_mm256_add_ps(_mm256_add_ps(valA, valB), valC), val025);

                __mmask16 up = _mm_cmpneq_epu8_mask(valXU, Zero128) | _mm_cmpneq_epu8_mask(_mm_loadu_epi8(YPtrU + j), Zero128);
                __mmask16 dn = _mm_cmpneq_epu8_mask(_mm_loadu_epi8(XPtrD + j), Zero128) | _mm_cmpneq_epu8_mask(_mm_loadu_epi8(YPtrD + j), Zero128);
                __mmask16 lf = _mm_cmpneq_epu8_mask(_mm_loadu_epi8(XPtr + j - 1), Zero128) | _mm_cmpneq_epu8_mask(valYL, Zero128);
                __mmask16 rt = _mm_cmpneq_epu8_mask(_mm_loadu_epi8(XPtr + j + 1), Zero128) | _mm_cmpneq_epu8_mask(_mm_loadu_epi8(YPtr + j + 1), Zero128);
                __mmask16 mask = up & dn & lf & rt;

                //_mm256_blend_ps(valZero, vallo, mask);

            }
        }
    }
}


void updateHeights(cv::Mat& Normals, cv::Mat& Z, int iterations) 
{
    for (int k = 0; k < iterations; k++) 
    {
        for (int i = 1; i < Normals.rows - 1; i++) 
        {
            for (int j = 1; j < Normals.cols - 1; j++) 
            {
                float zU = Z.at<float>(cv::Point(j, i - 1));
                float zD = Z.at<float>(cv::Point(j, i + 1));
                float zL = Z.at<float>(cv::Point(j - 1, i));
                float zR = Z.at<float>(cv::Point(j + 1, i));
                float nxC = Normals.at<cv::Vec3b>(cv::Point(j, i))[0];
                float nyC = Normals.at<cv::Vec3b>(cv::Point(j, i))[1];
                float nxU = Normals.at<cv::Vec3b>(cv::Point(j, i - 1))[0];
                float nyU = Normals.at<cv::Vec3b>(cv::Point(j, i - 1))[1];
                float nxD = Normals.at<cv::Vec3b>(cv::Point(j, i + 1))[0];
                float nyD = Normals.at<cv::Vec3b>(cv::Point(j, i + 1))[1];
                float nxL = Normals.at<cv::Vec3b>(cv::Point(j - 1, i))[0];
                float nyL = Normals.at<cv::Vec3b>(cv::Point(j - 1, i))[1];
                float nxR = Normals.at<cv::Vec3b>(cv::Point(j + 1, i))[0];
                float nyR = Normals.at<cv::Vec3b>(cv::Point(j + 1, i))[1];
                int up = nxU == 0 && nyU == 0 ? 0 : 1;
                int down = nxD == 0 && nyD == 0 ? 0 : 1;
                int left = nxL == 0 && nyL == 0 ? 0 : 1;
                int right = nxR == 0 && nyR == 0 ? 0 : 1;

                if (up > 0 && down > 0 && left > 0 && right > 0)
                {
                    Z.at<float>(cv::Point(j, i)) =
                        1.0f / 4.0f * (zD + zU + zR + zL + nxU - nxC + nyL - nyC);
                }
                else
                {
                    int stop = 0;
                }
            }
        }
    }
}

cv::Mat cvtFloatToGrayscale(cv::Mat F, int limit = 255) {
    double min, max;
    cv::minMaxIdx(F, &min, &max);
    cv::Mat adjMap;
    cv::convertScaleAbs(F, adjMap, limit / max);
    return adjMap;
}

cv::Mat localHeightfield(cv::Mat Normals) 
{
    const int pyramidLevels = 4;
    const int iterations = 200;

    std::vector<cv::Mat> pyrNormals;
    cv::Mat Normalmap = Normals.clone();
    pyrNormals.push_back(Normalmap);
    for (int i = 0; i < pyramidLevels; i++) 
    {
        cv::pyrDown(Normalmap, Normalmap);
        pyrNormals.push_back(Normalmap.clone());
    }


    cv::Mat Z(pyrNormals[pyramidLevels - 1].rows,pyrNormals[pyramidLevels - 1].cols, CV_32FC1, cv::Scalar::all(0));
    for (int i = pyramidLevels - 1; i > 0; i--) 
    {
        updateHeights(pyrNormals[i], Z, iterations);
        cv::pyrUp(Z, Z);
    }

    /* linear transformation of matrix values from [min,max] -> [a,b] */
    //double min, max;
    //cv::minMaxIdx(Z, &min, &max);
    //double a = 0.0, b = 150.0;

    //for (int i = 0; i < Normals.rows; i++) {
    //    for (int j = 0; j < Normals.cols; j++) {
    //        Z.at<float>(cv::Point(j, i)) =
    //            (float)a +
    //            (b - a) * ((Z.at<float>(cv::Point(j, i)) - min) / (max - min));
    //    }
    //}

    return Z;
}

//光源顺序 下 左 上 右
int mainlll(int argc, char* argv[]) 
{
    int screenWidth = 1440;
    int screenHeight = 900;
    int numPics = 4;


    __mmask16 a = 1;
    __mmask16 b = 4;

    __mmask16 c = a | b;


    std::vector<cv::Mat> camImages;

    for (int i = 0; i < numPics; i++) 
    {
        std::stringstream s;
        s << "images/image" << i << ".png";
        camImages.push_back(cv::imread(s.str(), IMREAD_GRAYSCALE));
    }

    cv::Mat S = computeNormals(camImages);
    cv::Mat XVal, YVal,  ZVal;
    NormalsSpeed(camImages, XVal, YVal, ZVal);



    cv::Mat Normalmap;
    cv::cvtColor(S, Normalmap, COLOR_BGR2RGB);
    cv::imshow("normalmap.png", Normalmap);
    /* compute depth map */
    cv::Mat Depth = localHeightfield(S);
    cv::imshow("Local Depthmap", cvtFloatToGrayscale(Depth));
    //exportMesh(Depth, S, camImages[0]);
    cv::waitKey(0);
    return 0;
}






void mainsda()
{
	Mat img1 = imread("nnn.jpg", 0);

	//Mat K,disX,disY;
	//getCannyKernel(K, 10);
	//Mat blur;
	//GaussianBlur(img1, blur, Size(0, 0), 20, 20);
	//tc.restart();
	//Filter2d(blur, K, &disX,&disY);
	//tc.stop();
	//std::cout << tc.elapsed_ms() << std::endl;
	//Mat d;
	//getCannyKernel(d, 10);
	//Mat one = Mat::ones(Size(1, 1), CV_16S);
	//Mat dx, dy;

	//tc.restart();
	//sepFilter2D(blur, dx, CV_16S, d, one);
	//sepFilter2D(blur, dy, CV_16S, one, d);
	//tc.stop();
	//std::cout << tc.elapsed_ms() << std::endl;

	Mat outimg1(img1.size(), CV_8UC1);
	Mat outimg2(img1.size(), CV_8UC1);

	Mat outimgX(img1.size(), CV_16SC1);
	Mat outimgY(img1.size(), CV_16SC1);

	std::vector<std::vector<cv::Point2f>> RList;
	tc.restart();
	
	EdgeSubpix(img1, 10, 35, RList, 4);
	tc.stop();
	std::cout << tc.elapsed_ms() << std::endl;

	RotatedRect rct;
	std::vector<cv::Point2f> pc;
	pc.push_back(Point2f(10,0));
	pc.push_back(Point2f(0, 10));
	pc.push_back(Point2f(-10, 0));
	pc.push_back(Point2f(0, -9));
	DirectFitEllipse(pc, &rct);

	outimg1.setTo(0);
	//cvtColor(img1, outimg1,COLOR_GRAY2BGR);

	for (auto item: RList)
	{
		//if (item.size()>100)
		//{
		//	DirectFitEllipse(item,&rct);

		//	ellipse(outimg1, rct, Scalar(0, 0, 255), 4);
		//}
		
		for (auto ite: item)
		{
			//outimg1.at<uchar>(ite.y, ite.x)[0] = 0;
			//outimg1.at<Vec3b>(ite.y, ite.x)[1] = 0;
			outimg1.at<uchar>(ite.y, ite.x) = 255;
		}
	}



	//HRegion reg;
	//thresholdRLE(img1, reg, 50, 150);

	SobelXY(img1, &outimgX);
	outimg1.setTo(0);
	ContourAnchor(outimgX, &outimg1, 25);


	Mat XBorder, YBorder, XYBorder;
	//Sobel(img1, XBorder, CV_16S, 1, 0, 3, 1.0, 0);
	//Sobel(img1, YBorder, CV_16S, 0, 1, 3, 1.0, 0);
	//convertScaleAbs(XBorder, XBorder);
	//convertScaleAbs(YBorder, YBorder);
	////XY方向上的因
	//addWeighted(XBorder, 0.5, YBorder, 0.5, 0, XYBorder);


	

	tc.restart();
	for (size_t i = 0; i < 1000; i++)
	{
		//SobelXY(img1, &outimgX, &outimgY);
		//Threshoud(img1, 80, 200, &outimg1);
		//threshold(img1, outimg1, 80, 200, 0);

		//Sobel(img1, XBorder, CV_16S, 1, 0, 3, 1.0, 0);
		//Sobel(img1, YBorder, CV_16S, 0, 1, 3, 1.0, 0);

		Threshoud(img1, 80, 200, &outimg1);
	}
	tc.stop();
	std::cout << tc.elapsed_ms() << std::endl;

	Sauvola(img1, &outimg1, 0.4, 300);

	FillHole(outimg1,true);

	//TotalVariationDenoising(img1, &outimg1);
	//DericheSpeed(img1, &outimg1, 0.07);



	//for (size_t i = 0; i < 100; i++)
	//{
	//	tc.restart();

	//	//Sauvola(img1, &outimg1, 0.4, 50);

	//	//TotalVariationDenoising(img1, &outimg1);
	//	DericheSpeed(img1, &outimg1, 0.07);

	//	tc.stop();
	//	std::cout << tc.elapsed_ms() << std::endl;
	//}


	system("pause");

}









